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Abstract 

We review some recent work in fast, efficient and accurate methods to 
compute viscosity solutions and non-viscosity solutions to static Hamilton- 
Jacobi equations which arise in optimal control, anisotropic front propagation, 
and multiple arrivals in wave propagation. For viscosity solutions, the class of 
algorithms are known as "Ordered Upwind Methods" , and rely on a systematic 
ordering inherent in the characteristic flow of information. For non-viscosity 
multiple arrivals, the techniques hinge on a static boundary value phase-space 
formulation which again can be solved through a systematic ordering. 
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1. Introduction 

This paper reviews recent work on algorithms for static Hamilton- Jacobi equa- 
tions of the form H(Du,x) — 0; the solution u depends on x S R n , and boundary 
conditions are supplied on a subset of R n . These equations arise in such areas as 
wave propagation, optimal control, anisotropic front propagation, medical imaging, 
optics, and robotic navigation. We develop algorithms to solve these equations re- 
markably quickly, with the same optimal efficiency as classic algorithms for shortest 
paths on discrete weighted networks, but extended to continuous Hamilton- Jacobi 
equations. 

The algorithms, which rely on a close examination of the flow of information 
inherent in static Hamilton- Jacobi equations, are robust, unconditionally stable 
without time step restriction, and efficient. They are "One-pass" schemes, in that 
the solution is computed at N grid points in 0(N log AT) steps. 
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1.1. Viscosity vs. non- viscosity solutions 

What is meant by a solution to H(Du, x) = 0? Viscosity solutions provide 
a unique and well-posed formulation which is linked to the unique viscosity limit 
of the associated smoothed equation; these are first arrivals in the propagation 
of information. Fig. shows an example from semiconductor manufacturing in 
which a beam whose strength is angle-dependent is used to anistropically etch away 
a metal surface. Fig. \%)p shows an optimal control problem to find the shortest exit 
path for a vehicle with position and direction-dependent speed. 




Ion etching in anisotropic front propagation 
Fig. Hi 




Optimal control 
Fig. 

Figure 1: Viscosity solutions to static HJ equations 



The above are viscosity solutions. However, there are many cases in which later 
arrivals, or "non-viscosity" solutions, are desirable. Fig. [2^ shows the propagation 
of a wave inwards from a square boundary; the evolving front passes through itself 
and later arrivals form cusps and swallowtails as they move; Fig. |2b shows multiple 
arrivals in geophysical wave propagation. 
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Fig. Eh 

Figure 2: Non- viscosity solutions 



Our goal is to create efficient algorithms which allow us to compute both types 
of solutions. In the case of viscosity solutions, algorithms are provided by the class 
of "Ordered Upwind Methods" developed by Sethian and Vladimirsky in |121 113| : 
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these methods work in physical space and construct the solution in a "One-pass" 
manner through a careful adherence to a causality inherent in the characteristic flow 
of the information. In the case of non-viscosity solution, algorithms are provided 
by the time-independent phase-space formulation developed by Fomcl and Sethian 
0, which relies on conversion of multiple arrivals into an Eulerian static boundary 
value problem, which can also be solved very efficiently in a "One-pass" manner 
which avoids all iteration through a careful ordering procedure. The remainder of 
this paper is devoted to describing these two classes of algorithms and providing a 
few computational results. 

2. Fast methods for viscosity solutions 

We first discuss "Ordered Upwind Methods" introduced in ^21 f° r computing 
viscosity solutions. 

2.1. Discrete control: Dijkstra's method 

Consider a discrete optimal trajectory problem on a network. Given a network 
and a cost associated with each node, the global optimal trajectory is the most 
efficient path from a starting point to some exit set in the domain. Dijkstra's 
classic algorithm 0] computes the minimal cost of reaching any node on a network 
in 0(N log N) operations. Since the cost can depend on both the particular node, 
and the particular link, Dijkstra's method applies to both isotropic and anisotropic 
control problems. The distinction is minor for discrete problems, but significant for 
continuous problems. Dijkstra's method is a "one-pass" algorithm; each point on 
the network is updated a constant number of times to produce the solution. This 
efficiency comes from a careful analysis of the direction of information propagation 
and stems from the optimality principle. 

We briefly summarize Dijsktra's method, since the flow logic will be important 
in explaining our Ordered Upwind Methods. For simplicity, imagine a rectangular 
grid of size h, where the cost dj > is given for passing through each grid point 
Xij = (ih,jh). Given a starting point, the minimal total cost Uij of arriving at 
the node be written in terms of the minimal total cost of arriving at its 

neighbors: 

Uij = mm(Ui-i t j,U i+ i t j,Ui t j-i,Ui t j + i) + Cy. (2.1) 

To find the minimal total cost, Dijkstra's method divides mesh points into 
three classes: Far (no information about the correct value of U is known), Accepted 
(the correct value ofU has been computed), and Considered (adjacent to Accepted). 
The algorithm proceeds by moving the smallest Considered value into the Accepted 
set, moving its Far neighbors into the Considered set, and recomputing all Con- 
sidered neighbors according to formula 12. II This algorithm has the computational 
complexity of 0(Nlog(N)); the factor of log(TV) reflects the necessity of maintain- 
ing a sorted list of the Considered values Ui to determine the next Accepted mesh 
point. Efficient implementation can be obtained using heap-sort data structures. 
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2.2. Continuous control: ordered upwind methods 

Consider now the problem of continuous optimal control; here, the goal is to 
find the optimal path from a starting position to an exit set. Dijkstra's method does 
not converge to the continuous solution as the mesh becomes finer and finer, since 
(see [lip it produces the solution to the partial differential equation max(|u s |, \u y \) = 
h*C, where h is the grid size. As h goes to zero, this does not converge to the solu- 
tion of the continuous Eikonal problem given by + Uy\ x ^ 2 = C. Thus, Dikstra's 
method cannot be used to obtain a solution to the continuous problem. 

2.2.1. Ordered upwind solvers for continuous isotropic control 

In the case of isotropic cost functions in which the cost depends only on po- 
sition and not on direction, two recent algorithms, first Tsitsiklis's Method ^J] 
and then Sethian's Fast Marching Method JUj have been introduced to solve the 
problems with the same computational complexity as Dijkstra's method. Both 
methods exploit information about the flow of information to obtain this efficiency; 
the causality allows one to build the solution in increasing order, which yields the 
Dijkstra-like nature of the solutions. Both algorithms result from a key feature 
of Eikonal equations, namely that their characteristic lines coincide with the gra- 
dient lines of the viscosity solution u(x); this allows the construction of one-pass 
algorithms. Tsitsiklis' algorithm evolved from studying isotropic min-time optimal 
trajectory problems, and involves solving a minimization problem to update the 
solution. Sethian's Fast Marching Method evolved from studying isotropic front 
propagation problems, and involves an upwind finite difference formulation to up- 
date the solution. Each method starts with a particular (and different) coupled 
discretization and each shows that the resulting system can be decoupled through 
a causality property. We refer the reader to these references for details on ordered 
upwind methods for Eikonal equations, as well as J3| for a detailed discussion about 
the similarities and differences between the two techniques. 

2.2.2. Ordered upwind solvers for continuous anisotropic general opti- 
mal control 

Consider now the full continuous optimal control problem, in which the cost 
function depends on both position and direction. In |12II13| . Sethian and Vladimirsky 
built and developed single-pass "Ordered Upwind Methods" for any continuous op- 
timal control problem. They showed how to to produce the solution Ui by recalcu- 
lating each Ui at most r times, where r depends only the equation and the mesh 
structure, but not upon the number of mesh points. 

Building one-pass Dijkstra-like methods for general optimal control is con- 
siderably more challenging than it is for the Eikonal case, since characteristics no 
longer coincide with gradient lines of the viscosity solution. Thus, characteristics 
and gradient lines may in fact lie in different simplexes. This is precisely why both 
Sethian's Fast Marching Method and Tsitsiklis' Algorithm cannot be directly ap- 
plied in the anisotropic (non-Eikonal) case: it is no longer possible to de-couple the 
system by computing/accepting the mesh points in the ascending order. 
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The key idea introduced in ^1 EH is to use the location anisotropy of the 
cost function to limit of the number of points on the accepted front that must be 
examined in the update of each Considered point. Consider the anisotropic min- 
time optimal trajectory problems, in which the speed of motion depends not only 
on position but also on direction. The value function u for such problems is the 
viscosity solution of the static Hamilton- Jacobi-Bellman equation 



In this formulation, a is the unit vector determining the direction of motion, f(a, x) 
is the speed of motion in the direction a starting from the point x G O, and q{x) 
is the time-penalty for exiting the domain at the point x G dfl. The maximizer a 
corresponds to the characteristic direction for the point x. If / does not depend on 
a, Eqn. 12.21 reduces to the Eikonal equation, see pp. 

Now, define the anisotropy ratio F1/F2, where < F\ < f(a,x) < F2 < 00. 
In |13|. two key lemmas were proved: 

• Lemma 1. Consider the characteristic passing through x G Q and level curve 
u(x) = C, where q m ax < C < u(x). The characteristic intersects that level set 
at some point x. If x is distance d away from the level set then \\x — x\\ < 

• Lemma 2. Consider an unstructured mesh X of diameter h on 51. Consider 
a simple closed curve T lying inside Q with the property that for any point x 
on F, there exists a mesh point y inside T such that \\x — y\\ < h. Suppose the 
mesh point Xi has the smallest value u(x~i) of all of the mesh points inside the 
curve. If the characteristic passing through of, intersects that curve at some 
point £i then \\£i — Xi\\ < h^r. 

Thus, one may use the anisotropy ratio to exclude a large fraction of points on 
the Accepted Front in the update of any Considered Point; the size of this excluded 
subset depends on the anisotropy ratio. Building on these results, a fast, Dijkstra- 
like method was constructed. As before, three of mesh points classes are used. The 
Accepted Front is defined as a set of Accepted mesh points, which are adjacent to 
some not-yet-accepted mesh points. Define the set AF of the line segments XjXk, 
where Xj and Xk are adjacent mesh points on the AcceptedFront, such that there 
exists a Considered mesh point Xi adjacent to both Xj and Xk- For each Considered 
mesh point Xi one defines the part of AF "relevant to Xi" : 



We will further assume that some consistent upwinding update formula is available: 
if the characteristic for xi lies in the simplex XiXjXk then Ui = K(Uj, Uk,Xi, xj, Xk)- 
For the sake of notational simplicity we will refer to this value as Kj^. 

1. Start with all mesh points in Far {Ui = 00). 

2. Move the boundary mesh points (xi G Sil) to Accepted (Ui = q(xi)). 

3. Move all the mesh points x% adjacent to the boundary into Considered and 
evaluate the tentative value of Ui = min( x . iX k)N Ffe) Kj,k- 



max aeSl {(V«(x) • (-o))/(o, x)} = 1, x £ CI, 
u(x) = q(x), x G 30. 



(2.2) 




(xj,Xk) G AF \3x on (xj,Xk) s.t. ||x - if|| < h— 
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4. Find the mesh point x r with the smallest value of U among all the Considered. 

5. Move x r to Accepted and update the Accepted Front. 

6. Move the Far mesh points adjacent to x r into Considered. 

7. Recompute the value for all the Considered Xi within the distance hjjr from 
x r . If less than the previous tentative value for Xi then update £/;. 

8. If Considered is not empty then go to 4). 



2.2.3. Analysis and results 

This is a "single-pass" algorithm since the maximum number of times each 
mesh point can be re-evaluated is bounded by the number of mesh points in the 
neighborhood of that point; the method formally has the computational complexity 
ofO((^) 2 Mlog(M)). Conver gence of the method to the viscosity solution is proved 
in [T3|, and depends on the upwinding update formula Ui — K(Uj, Uk, Xi, Xj , Xk). 

As an example, taken from we compute the geodesic distance on the 
manifold g(x,y) = .9sin(27rx) sin(27n/) from the origin. This can be shown to be 
equivalent to solving the static Hamilton- Jacobi equation 



\\Vu(x)\\F[x, 



Vu(x) 
\Vu(x)\ 



= 1, 



(2.3) 



-+<?y cos 2 (w) + (? 2 sin 2 (oj)— g x g v sin(2oj) 



with speed function F given by F(x,y,u>) — y 1+g2+g2 

where ui is the angle between Vu(x, y) and the positive direction of the x-axis. The 
anisotropy is substantial, since the dependence of F upon uj can be pronounced 
when Vg is relatively large. Equidistant contours are shown on the left in Figure 




i Want Arrival Path 




SO 100 150 200 250 300 3 50 



Source 

Figure 3: Left: Anistropic front propagation Right: Arrival paths 
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3. Fast methods for multiple arrivals 

3.1. Computing multiple arrivals 

We now consider the problem of multiple arrivals. As an example, consider 
the two-dimensional Eikonal equation 

\Vu\F(x,y) = (3.1) 

with F(x,y) given. We imagine a computational domain and a source point, as 
shown on the right in Figure |31 Suppose the goal is to determine the arrival time 
and path to each point in the interior from the source point. Here, we are interested 
not only in the first arrival, but all later arrivals as well. 

One popular approach to computing multiple arrivals is to work in phase 
space, in which the dimensionality of the problem is increased from physical space to 
include the derivative of the solution as well. There are two approaches to computing 
these multiple arrivals through a phase space formulation. One is a Lagrangian (ray 
tracing) approach, in which the phase space characteristic equations are integrated, 
often from a source point, resulting in a Lagrangian structure which fans out over 
the domain. Difficulties can occur in either in low ray density zones where there 
are very few rays or near caustics where rays cross. The other is an Eulerian 
description of the problem, in either the physical domain or phase space. In recent 
years, this has led to many fascinating and clever Eulerian PDE-based approaches 
to computing multiple arrivals, see, for example, ^3 El El 03 G] • We note that the 
regularity of the phase space has been utilized previously in theoretical studies on 
the asymptotic wave propagation [H] . The above phase space approaches to solving 
for multiple arrivals have two characteristics in common: 

• A phase space formulation increases the dimensionality of the problem. In two 
physical dimensions, the phase space formulation requires three dimensions; in 
three physical dimensions, the phase space formulation is in five dimensions. 

• Given particular sources, the problem is solved with those source location(s) 
as initial data. Different sources requires re-solving the entire problem. 

• The problem is cast as an initial value partial differential equation, and is 
evolved in time. Time step considerations in regions of high velocity play a 
role in the stability of the underlying scheme. 

3.2. A boundary value formulation 

Fomel and Sethian [7J take a different approach. A set of time-independent 
"Escape Equations" are derived, each of which is an Eulerian boundary value partial 
differential equation in phase space. Together, they give the exit time, location and 
derivative of all possible trajectories starting from all possible interior points. Thus, 
the particular choice of sources is reduced to post-processing. The computational 
speed depends on whether one wants to obtain results for all possible boundary 
conditions, or in fact only for a particular subset of possibilities. 
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3.2.1. Liouville formulation 

Briefly (see for details) begin with the static Hamilton- Jacobi equation 

ff(i,Vii) = 0, (3.2) 

and write the well-known characteristic equations in phase space (x,p), where p 
corresponds to Vit (see, for example, 6 ). The characteristics must obey 

= V pJ ff; ^ = -V X H. (3.3) 

da da 

Differentiating the function u(x(a)), we obtain an additional equation for trans- 
porting the function u along the characteristics: 

du _ dx _ TT . 

— = Vu ■ — = p ■ V P H. (3.4) 
da da 

Eqns. 13.313.41 can be initialized at a = 0: x(0) = Xo, p(0) = po, u(0) = 0. 

One can now convert the phase space approach into a set of Liouville equations. 
To simplify notation, we denote the phase-space vector (x,p), by y, the right-hand 
side of system given in Eqn. 13.31 by vector function R(y), and the right-hand side 
of Eqn. 13.41 by the function r(y). In this notation, the Hamilton-Jacobi system is 

d -^ = R(y); 9 -^ = r(y), (3.5) 

oa da 
and is initialized at a — as y = yo and u = 0. This system satisfies 

— — — = V 2/i?(y ) , (3.6) 
and the transported function u satisfies the analogous equation 

^^ = Vo.^o) + r(y ), (3.7) 
where Vq denotes the gradient with respect to yo. These are the Liouville equations. 



3.2.2. Formulation of escape equations 

The key idea in |Jj is as follows. Assume a closed boundary &D in the y 
space that is crossed by every characteristic trajectory originating in yo £ T>. This 
defines for every yo the function a = a(yo) of the first crossing of the corresponding 
characteristic with dT>. Now introduce a differentiable function T(y) that identifies 
the boundary, that is, T(y) = 0. In particular, we then have that T {y{yo, &{yo)) = 0. 
One can then differentiate with respect to the initial condition y to obtain an escape 
equation for the parameter a. Similarly, one can derive escape equations for the 
position and value, yielding the full set of 



Escape Equations 1 + Voa ■ R(yo) = 

S7 yR(yo) = o 

V u- R(yo)+r{y ) = 
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3.3. Fast solution of escape equations 

Summarizing, rather than compute in physical space, we derive boundary value 
Escape equations in phase space y = (x,p). All time step considerations are avoided, 
and one can compute all the arrivals from all possible sources simultaneously. This 
Eulerian formulation means that the entire domain is covered, even quiet slow zones. 

Finally, and most importantly, a constructive, "One-pass" algorithm, similar 
to the one presented for viscosity solutions, can be designed. Exit time, position, and 
derivative at the boundary form boundary conditions. We can then systematically 
march the solution inwards in phase space from the boundary, constructing the 
solution through an ordering sequence based on the characteristics that ensures 
computational phase space mesh points need not be revisited more than once. 

Consider a square boundary as an example, and suppose we wish to find the 
time u(x, z, 9) at which a ray leaving the initial point (x, z) inside the square, initially 
moving in direction 9, hits the boundary. We assume that the slowness field n(x, z) 
is given. First, note that the set u(x, z, 9) — T, drawn in x, z, 9 space, gives the 
set of all initial positions and directions which reach the boundary of the square at 
time T. By the uniqueness of characteristics, the set of all points parameterized 
by T and given by U(T) = {x, z, 9 \ u(x, z, 9) — T} sweep out the solution space. 
Figure^ shows the solution surfaces u(x, z, 9) for the collapsing square. 

Details on the exact algorithm are given in [J]. As demonstration (see jjj), in 
FigureEb, the top pair shows all the arrivals starting from a source at the center of 
the top wall, together with the slowness field on the right (darker is slower). The 
bottom pair shows the first arrival and on the amplitude of the displayed arrival 
(the lighter the tone, the more amplitude). 
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